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ABSTRACT 


Breather solitons in a one-dimensional lattice of coupled nonlinear oscillators are 
numerically investigated. These are localized nonpropagating steady states that exist at 
frequencies either below the linear cutoff frequency (corresponding to the extended mode in 
which all the oscillators are in-phase) or above the upper linear cutoff frequency (corresponding 
to the extended mode in which each oscillator is 180° out-of-phase with its immediate neighbors). 
The lattice is damped and parametrically driven. A nonlinear Schrodinger theory, which assumes 
a modulational amplitude that is weakly nonlinear and slowly varying in space, is compared to 
the numerical data. The error is roughly 5% at low amplitudes and 20% at high amplitudes. 
The regions in the drive parameter plane (amplitude vs frequency) where the breathers exist are 
numerically determined and compared to theory. A substantial discrepancy occurs at lower drive 
amplitudes where the theory predicts that the lower cutoff breather should exist, but where an 
instability is observed. Also in contrast to the theory, the region of the upper cutoff breather has 
relatively large areas in which quasiperiodicity occurs or the motion decays to rest. 
Quasiperiodicity is also observed in the lower cutoff breather. Finally, instead of a global 
parametric drive, an end drive is investigated. It is found that, for drive frequencies outside the 
linear propagation band, there is an amplitude threshold for the periodic ejection or "shedding" 
of propagating breather solitons. The quasiperiodicity that occurs for a global parametric drive 


may be a consequence of soliton shedding. 
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I, INTRODUCTION 


A. BACKGROUND 


A soliton is a localized wave whose shape or envelope is constant as a result of 





nonlinearities and dispersion. Furthermore, it collides elastically with other waves. Solitons 
were first observed as canal surface waves by John Scott Russell in 1834 (Dodd et al. 1978). 
Other systems are now known to support solitons of various types. Research in fiber optic 
solitons is currently very active because of their potential use in communications and optical 
switching (Mollenauer et al. 1991). 

In this thesis, a one-dimensional lattice of coupled nonlinear oscillators is numerically 
investigated. This system is known to possess kink solitons (Galvin 1990, Denardo et al. 1991). 
It will be shown that breather solitons can also exist. These are nonpropagating steady states in 
which most of the oscillators are approximately at rest while those in a localized region have 
substantial amplitude. Such a state is clearly not a linear mode of the system. Breathers are self- 
trapped states that occur at frequencies outside the linear propagation band. Two types exist: 
lower cutoff breathers, in which the oscillators are in-phase, and upper cutoff breathers, in which 
the oscillators are 180° out-of-phase. Lower and upper cutoff solitons can be considered as 
finite-amplitude modulations of the uniform lower and upper cutoff modes, respectively. 

In the absence of drive and dissipation, it is known that weakly nonlinear breathers 
described by the ¢* theory are unstable, although the decay rate is extremely slow (Segur and 
Kruskal 1987). It is also known that instabilities can occur as a result of discreteness in systems 
that are undriven and undamped (Goedde 1990). In this thesis, dissipation and global drive are 


employed, which lead to stable breathers. 


The simplicity of the model equation suggests that breathers can occur in a variety of 
systems. Indeed, the first observation of breathers was as cross surface waves on a long channel 
of liquid (Wu et al. 1984). Lower cutoff breathers have been observed in a pendulum lattice 
(Denardo 1990). Recently, upper cutoff breathers have been observed in a magnetically coupled 
pendulum lattice and confirmed numerically with an equation that models this system (Atchley 
1991). 

It is also shown that, for a local pure-frequency drive confined to one end site of a lattice, 
it is possible to shed propagating solitons. The shedding occurs at a frequency that is 
quasiperiodic relative to the drive frequency. This phenomenon has previously been observed 
in only one system, surface waves just below the second cutoff mode in a large tank (Kit et al. 
1987, Shemer 1990). The existence of the soliton shedding in a simple lattice suggests that it is 
a general phenomenon which can occur in a variety of systems, and that a simple fundamental 


explanation should exist, although no such explanation is known at present. 


B. EQUATION OF MOTION 

To motivate the equation of motion that will be investigated, a model system is considered. 
The system is a uniform lattice of linearly coupled simple pendulums that oscillate transverse to 
the lattice. The torque on the nth pendulum is: 


T= he (0 200 ee Oem Te Bach: 


where p is the torsional constant that characterizes the coupling, m is the pendulum mass, and 


| is the pendulum length. For weakly nonlinear motion, the approximation 


sin(@) =6 - 2 63 TBE 


sin(§) =6 - 2 9° Ten 


can be made. Linear damping is included by adding to I.B.1 a torque proportional to the velocity 
6’, (with a negative coefficient), where the prime denotes differentiation with respect to time. 
To sustain localized structures in the lattice, a global parametric drive that results from vertically 
oscillating the lattice, is utilized. By the Equivalence Principle, in the reference frame of the 


lattice the effective acceleration due to gravity is: 


Jesp FJ + A cos(2ut), 1G 03} G: 


where a is the acceleration amplitude of the drive and 2w is the drive frequency. 
Combining the above effects, and using Newton’s Second Law, the equation of motion of 


for the lattice can be expressed as: 


67,-c?(6,,, - 20, + 6,.,) +B0’,+[w ,2+2ncos(2wt)]6, = «8,7, I.B.4 
where c* is a measure of the coupling between lattice sites and w, is the cutoff frequency of a 
pendulum. For simplicity, the term proportional to the product of the drive and nonlinearity has 
been dropped. It can be shown that this has no essential effect upon the localized structures. The 
nonlinear coefficient @ equals w,’/6 if (I1.B.4) is to approximate a pendulum lattice. The model 
is generalized by allowing the nonlinear coefficient to be positive or negative. The lattice is 
driven parametrically with an amplitude of 27. The fundamental response frequency w is half the 
drive frequency. The lower cutoff mode of the lattice is characterized by uniform motion in the 


lattice (Fig. I.B.1), and has linear frequency w,. The upper cutoff mode has a 180° phase 


difference between adjacent lattice sites (Fig. I.B.2), and has linear frequency (w,? + 4c?)'*. For 


a softening (a > Q) lattice, the upper cutoff mode is stable. The lower cutoff mode is subject 
to the Benjamin-Feir instability (Denardo 1990), and the motion evolves into one or more 
breather solitons (Fig. I.B.3). For a hardening (@ < 0) lattice, the roles of the cutoff modes are 


reversed and an upper cutoff breather evolves (Fig. 1.B.4). 


HF 


Figure I.B.1 Lower Cutoff Mode 


HTT 


Figure I.B.2 Upper Cutoff Mode 





Figure I.B.3 Lower Cutoff Breather 





Figure I.B.4 Upper Cutoff Breather 


Il. THEORY 


A. LOWER CUTOFF BREATHER 
For modulations of the lower cutoff mode, the displacement (6,) in I.B.4 can be written 
approximately as: 
8 -A(n, t) e7°°+ complex conjugate, Tear 


where A(x,t) is a complex differentiable function of its arguments. The lattice spacing is assumed 
to be unity. A displacement (6,) containing a third harmonic term produced approximately the 
same results aS II.A.1 and was consequently discarded. Substituting II.A.1 into 1.B.4 and 


requiring that A be slowly varying and weakly nonlinear, one obtains the modulational equation 


2iwA - c? A, + (w,2 - wo? + iwB)A+ 1 A* = BA lLAP A, TI fee 


which is a nonlinear Schrédinger equation. If the modulation is assumed to be nonpropagating, 


a Stationary solution 


A (x, £) = ANG) "en Tei ane 
exists where A and 6 are real. Substituting II.A.3 into II.A.2 gives: 
- c? A” + (w,? - w?)A - 3a A? + iwPA = -nAe??, IT Japa 


Equating the imaginary parts of both sides yields an expression for 6: 


sin(26) = oP ies 


The real part of I].A.4 yields: 


- iA sa 0, ig aS 


where BP - O.? - o2 + Vn? - wf. 


There are actually two possibilities because of the radical. However, for a softening lattice the 
negative branch leads to an unstable solution (Denardo 1990). 
Equation II.A.6 can be integrated to yield an elliptic function. A localized solution to 


I].A.6 is of the form A = a sech(bx), provided 





Then A~-,|“H sech sei AG 


3a 








, Bese (x, - i) 
c? 


Substituting the solution for A into II.A.1, the expression for the displacement is: 


6. -2,)24 sech 2 Ca. ees (Geet 10!) JOE ae 
G 


This solution is valid for p > QO and a > QO (softening). 

The physical reasoning for the existence of the breather is derived from the curvature of 
the modulation envelope. For the softening lower cutoff breather, negative curvature in the body 
of the envelope produces a restoring force (Fig. II.A.1). The frequency in the body is below 
cutoff due to the nonlinearity which decreases the frequency more than the curvature increases 
it. In the tail of the envelope the curvature is positive, thus an anti-restoring force. The 


frequency is below cutoff and the breather is evanescent in the tail. At the inflection points, the 


nonlinearity is the sole reason the frequency is below cutoff. The breather is thus a self-trapped 

state due to the combination of the nonlinearity and curvature of the modulation. As a result, 

Steady state motion is possible. This motion has been observed in the simulated lattice. 
Concluding that steady state motion is possible, one should be able to determine the values 


of 7 and w where a stable solution exists. Requiring p» from II.A.6 to be positive, 

n?> (@,? - wo? )* + wo? B?. “II.A.9 
If w=w),, this yields the Matthieu hyperbola. The sech solution (Eqn. II.A.8) is unstable inside 
the hyperbola because of growth in the wings. There are two conditions for the solution to exist: 


n > wand p > 0. The first condition is 7 > w,6 if w=w,. The second condition is related 


to If.A.9, which can be rewritten as 


J 7 = @2) Paleo ie el, Ape 


If w > w,, then p > 0. However, if w < w, the requirement for np > 0 is 7 > wf. The sech 


solution decays to rest below the boundary 7 > wf. 


NEGATIVE 


POSITIVE 


POSITIVE 





Figure II.A.1 Lower Cutoff Breather Envelope Curvatures 


B. UPPER CUTOFF BREATHER 


The displacement (0,) for modulations of the upper cutoff mode can be written as 


6 = (-1) 7 A(n, t) e7°*+ complex conjugate Emre se I 


n 
Where A(n,t) is a complex differentiable function of its arguments. Again, the higher order terms 
have been neglected, and the lattice spacing is assumed to be unity. The equation that describes 


the weakly nonlinear complex amplitude is again a nonlinear Schrédinger equation: 


2iwA, + c? A, + (w,? - wo? + iwPB)A+ NA - 30/1 AP A, IT Bee 


Equation II.B.2 is identical to I].A.2 with the substitutions w,2-w,? and c*>-c”, where w,? = @9° 
+ 4c’, which is the square of the linear frequency of the upper cutoff mode. Assuming a 


nonpropagating modulation, a stationary solution for A Is: 


A.(e,t) =8A Gone Tee ef 


where A(x) and 6 are real. Substituting II.B.3 into II.B.2 gives: 
c? A” + (w,? - w?)A - 3a A? + iwPA = -nAe~2?, II.B.4 
Setting the imaginary parts of both sides equal and solving for 6, the phase relation is again: 


sin(26) = <r. TEASAS 


The real part of II.B.4 yields: 


-c? AYWY+vA+3a4A7=0 II.B.6 


Where v= 0? - w,? + ¥n* - w*f? 


10 


which is similar to I1.A.6. Again, the selected sign of the radical leads to a stable solution. A 


trial solution of the form A = a sech(bx) used in II.B.4 yields: 


Vv 
\3 (os) 


Substituting the solution for A into 11.B.1, the displacement (©,) for the upper cutoff mode is: 


0.7 (ayn | 2 sech \ CPO 
oi ce 


which is valid for y>0O and a< 0 (hardening). 





sech : 1 Bad, 








cos(wt + 6), II.B.8 








The physical mechanism for the existence of the upper cutoff breather is essentially 
identical to that of the lower cutoff breather. The negative curvature in the body of the hardening 
upper cutoff breather now produces an anti-restoring force (Fig. II.B.1) (Denardo 1990), 
however, the nonlinearity is stronger than the curvature and the frequency in the body is above 
the linear upper cutoff frequency. The positive curvature in the tails of the modulation creates 
a restoring force and the frequency is again above the cutoff frequency. Similar to before, the 
frequency is above the cutoff frequency at the inflection points as a sole result of the nonlinearity. 
Thus, the upper cutoff breather is also a self-trapped state due to the combination of curvature 
and nonlinearity. 

As for the lower cutoff case, a region of stable, steady state motion should also exist for 


the upper cutoff breather. Constraining v from II.B.6 to be positive, 


n@ > ( w,? - w? )? + w? B?. II.B.9 


11 


Again, if w~w, this yields the Matthieu hyperbola. The wings of the sech solution are unstable 
inside the hyperbola. There are two conditions for the solution exist: 7 > w@ and vy > 0. The 
first condition is 7 > w,8 if w=w,. The second condition is related to II.A.20 which can be 


rewritten as 


Vn? - wo? B? > Iw,? - oA. II.B.10 
If w > w,, then p > 0. However, if w < w, the requirement for p > 0 is 7 > wf. The sech 


solution decays to rest below the boundary 7 > w. 


2 


POSITIVE POSITIVE 





Figure II.B.1 Upper Cutoff Breather Envelope Curvatures 
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C. END DRIVEN LATTICES 

Nonlinear lattices driven at one end and at a frequency outside the linear propagation band 
have been observed to periodically eject solitons (Sect.III.D). This is in strong contrast to the 
linear theory which predicts only a steady state evanescent wave, as will now be shown. The 
equation of motion for a lattice of linearly coupled nonlinear oscillators (1.B.4) can be modified 
to replace the parametric drive with a direct drive. Further restricting the drive to the end site 


alone and removing the periodic boundary condition yields: 


n=0: 05-c?(0, - 0, )+PO,+w,2 8,+ncos(wt) = «6,2 Ll. Cy 


nO} e-c2(0., - 6.) +B0’ +o .? 6, - «6,>. 


n+1 
This produces an end driven lattice. The lattice is still linearly coupled and nonlinear, but is now 
more suitable for the observation of propagating solitons which have a spatially varying phase 
(Denardo 1990). Focusing on the linear motion in the continuum limit, equations II.C.2 have 
the form: 

be 0) - Ver - C7 Vi, + BY, + ©? y= 0 Lees 


ait (6)  - Ver *—— Vy +t OL? V+ HN Cos(Wt). 


Consider a solution of the form: 


y= A eter ex | Conplemcony uGaner Gy ee Geer 


with a complex wavenumber x=k+ia. Substituting II.C.4 into II.C.3 yields: 


-w* - c* (ik - a)* - ioB + w,? - 0. LERGaS 


Separating the real and imaginary parts of II.C.5 gives: 


Im: ee LOB 


2c? 


The case of interest is w@<,. The solution for B=0 is: 


k=0 and a = ,| —“——_ SOS ais 


Physically, the wave is purely exponential in space; the sign of the displacement is the same for 
every lattice site. If 640 then k can be eliminated using the real and imaginary parts of II.C.5 
to obtain: 
wD 2 — @2 - 
i op = Q. WS PS Gary 
Ze ZC" 


Solving for a and choosing the positive root , since @ was assumed to be real, gives: 


ale 
26. 





jo 2 - 2 + f(o,? - 07)? + (wB)?]. TievGae 


2 = 
a fe) Qo) 


The solution for k is then: 





oe oF i. II.c.9 


G@ (rs G)? 
a =-,{|—————-_ and k- OB __i__. Tt. G10 


15 


The damping gives rise to a spatial phase (k # 0). Specializing to strong damping, w8 > w,’-w* 


Or near-propagation, yields: 


ao k- | ob II.C.11 
2c? 


This is the largest that k/a can be (unity). In general, O<k/a<1 (Fig. II.C.1). 
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Figure II 
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Ill. NUMERICAL SIMULATION 


A. IMPLEMENTATION 

The motion of the lattice of linearly coupled nonlinear oscillators was solved with the 
Euler-Cromer method. This implementation utilized acceleration as the basic numerical iteration. 
The position calculation of each lattice site after a given time step was divided into three sections: 

1. The calculation of acceleration for each lattice site due to the position of it and the 

two adjacent sites. 

2. The calculation of the new velocity for each site by multiplying the acceleration by the time 

step, and adding the old velocity. 

3. The calculation of the new position for each lattice site using the new velocity for each site. 
Periodic boundary conditions were imposed, thus making the lattice a ring lattice with any 
curvature ignored. The Euler-Cromer method converged, with an increasing number of time 
Steps, to a "true" lattice amplitude (Fig. III.A.1). A time step that gave an amplitude within .1% 
of the "true" response amplitude allowed the program to run at a fairly high rate of speed with 
minimal loss of accuracy. 

The program is highly interactive and gives a "real time" picture of the lattice motion. The 
numerical simulation imposes some interactive restrictions. The elapsed time in the model frame 
must be reset periodically (preferably every period of the drive) to avoid strong transients when 
changing the drive frequency. The problem stems from the parametric drive term (cos(2wt)) 
Since the time increment is based on the response period. Thus a change in frequency in the 
middle of the cycle produces a discontinuity in the size of the time increment. The solution is 


to employ an integer counter to reset the elapsed time after one period of the drive and to require 


18 


any frequency changes to occur when the time is reset. The startup of a modulation envelope or 
"profile" can be accomplished in two ways: restarting a saved profile, or starting a profile created 
from theory. A profile can be saved while the program is running, thus eliminating the need to 
duplicate steps required to reach a particular set of parameters. All necessary parameters such 
as drive amplitude and frequency in addition to the position and velocity for each lattice site are 
written to a file of the researcher’s choice. The modulation must meet certain stability criterion 
to be recorded and is "captured" at the upper turning point (response) of a selected site. A 
previously saved profile can be recalled using only the filename. Since the response and drive 
have a phase difference, the program restarts the profile with a time phase correction to minimize 
startup transients. A startup file could also be generated using one of the programs that calculate 
a profile from theory using keyboard input of drive, coupling and other parameters. The 
equations of motion are in a subroutine, allowing them to be altered or replaced quickly. A 
different method of calculating the lattice site positions would probably not cause any difficulty 
as long as the new coordinates were available at the same place in the program as before. A 
complete list and more detailed description of the program options as well as a program code 


listing are included in the Appendix. 
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Figure III.A.1 Response Accuracy Plot 
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B. LOWER CUTOFF BREATHERS 

In order to investigate the stability of lower cutoff breathers, a drive plane, consisting of 
drive amplitudes and frequencies where the particular structure exists, was obtained by numerical 
solution of the equation of motion. The structure chosen for the lower cutoff drive plane was a 
symmetric breather (Fig. III.B.1) in which the maximum amplitude is centered on a lattice site 
(on-site). The nonlinear parameter (a) was 1 (lower cutoff), the coupling (c? or y) was .1 
(weak), dissipation (@) was .03 (small) and the lattice consisted of 50 sites long. The points in 
the drive plane (Fig. III.B.2) were almost exclusively obtained by incrementing the drive 
amplitude by .1% (or less) and recording the amplitude at which transitions occurred. 
Additionally, to probe stability, a <+3% (of the maximum response amplitude in the lattice) 
random kick was applied to the lattice after each increment. 

The upper boundary from w=.76 to .89 was characterized by the lattice "going over the 
top’ after perturbation (random kick). Because, the potential well for the lower cutoff mode 
(Fig. III.B.3) is not infinitely high, the lattice is able to escape the well (go over the top) and 
become unstable. There were no visual indications in the lattice of an impending instability 
before the boundary was reached. The right upper boundary, from w=.89 to .94, was marked 
by the generation of another symmetric breather (out of phase with the first) on the opposite side 
of the lattice ring from the initial breather. The lattice then went over the top from a random site 
fairly quickly after the generated breather reached an amplitude comparable to the original 
breather. Further tests with a lattice of 100 sites indicated that the generated breather is 
approximately 27 sites away from the center of the original breather and will arise on both sides 
of the original if the sites are available. The breather structure in the extreme lower right corner 


is quasiperiodic until the boundaries where it goes over the top without the generation of another 
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breather. The lower boundary from w=.85 to .947 is preceded by a quasiperiodic region. 
Ultimately the lattice either goes over the top or decays to rest after very strong transients. 
Careful observations are suggestive of quasiperiodicity resulting from the energy loss due to 
soliton shedding. The left lower boundary from w=.76 to .85 was only preceded by the 
quasiperiodic region at the points indicated but the strong transients were consistently observed 
before the lattice went over the top. A comparison of the stable region predicted by theory for 
the lower cutoff symmetric breather shows good agreement for the upper right boundary but a 
discrepancy on the lower boundary (Fig. III.B.4). Some destabilizing factor has prevented the 
lower boundary from reaching the theoretical limit. Experimentation with a half-site centered 
symmetric breather (Fig. III.B.5) revealed a link between the strong transients observed in the 
on-site breather drive plane but did not extend the lower boundary to the theoretical limit. The 
strong transients occur at the same drive amplitudes that the half-site breather shows quasiperiodic 
behavior. Increasing the drive amplitude results in a transition from half-site to on-site at 
approximately the start of the quasiperiodic line in the on-site drive plane. The half-site breather 
is stable as drive amplitude is decreased, but it decays to rest at approximately 7=0.9 for w=0.9. 
The transition from steady state motion to rest is interesting in that it shifts back to an on-site 
breather before it decays to rest. 

Comparison of the theoretical profile for the parameters used in the drive plane was done 
for the lowest and highest amplitude saved profiles as well as a half-site profile. The lowest 
amplitude comparison was sharper and lower in amplitude than the predicted profile (Fig. 
III.B.6). The highest amplitude comparison was also sharper and even lower in comparable 
amplitude (Fig. III.B.7). A comparison of the half-site profile shows much closer agreement with 
theory although the data still exhibits a greater amplitude (Fig. IJI.B.8). The amplitude of the 


saved profile from the program was not exactly the maximum amplitude of that particular state 
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did not precisely coincide with the peak response amplitude. This was corrected by using the 
maximum coordinate plus the coordinates for the preceding and following time steps to find the 
maximum amplitude by parabolic curve fitting. By recording the time step, the data could be 


scaled according to the time difference between the turning point and actual capture time. 
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Figure III.B.5 Symmetric Half-Site Lower Cutoff Breather 
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Figure III.B.7 Lower Cutoff Breather Theory Comparison 
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Figure III.B.8 Lower Cutoff Breather Theory Comparison 
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C. UPPER CUTOFF BREATHERS 

The structure chosen for the upper cutoff drive plane was the symmetric on-site breather 
(Fig. I1I.C.1). The parameters for the drive plane were a=-1 (hardening), y=.1 (weak 
coupling), 6=.03 (small) and a lattice length of 50 sites. The points in the drive plane (Fig. 
III.C.2) were almost exclusively obtained by incrementing the drive amplitude by .001 and 
recording the amplitude at which instabilities were observed. After each increment, a random 
kick of <+3% (of the maximum lattice response amplitude) was applied to the lattice. 

The upper left boundary marks the sharp transition from a breather structure to a complex 
upper cutoff mode structure. The transition time was fairly short and was characterized by the 
growth of one section of the lattice to an amplitude comparable to that of the original breather. 
The growing section was in a pure upper cutoff mode and spread to include the rest of the lattice 
after reaching full amplitude. The area inside the circles (Q) was where the lattice motion was 
quasiperiodic. The boundary was difficult to ascertain for omega 1.25 and higher since local 
regions appeared to exhibit quasiperiodicity but in fact merely had long settling times. Regions 
marked H denote areas of no lattice motion. Inside the region and below the lower boundary, 
the lattice decayed to rest with minimal transients. Region QS was characterized by a 
quasiperiodic shedding breather. The quasiperiodic breather was usually observed near the 
boundaries of QC and was similar to Fig. III.C.1. The center site of the breather was 
quasiperiodic and small solitons were ejected from the center in both directions simultaneously. 
The quasiperiodic breather spontaneously transitioned to a complex cyclic state at the lower 
boundary of QC (asterisks). The lattice motion in the small asterisked area (QC) was of two 
general types, a complex cyclic state or an anti-symmetric half-site breather. The cyclic state that 


evolved from the quasiperiodic breather appeared to be chaotic and then settled into multiple on- 
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site breathers. The motion shifted back and forth between the two states with a fairly long cycle. 
If the lattice was in this particular mode when the boundary to H (right side of QC) was reached, 
the cycle was broken and the lattice settled into a multiple breather state. The multiple breathers 
merged two by two until a single on-site breather was left. The single breather then decayed to 
rest. The lattice also spontaneously shifted from the complex cycle to a half-site antisymmetric 
breather (Fig. III.C.3). A half-site breather consistently evolved from an quasiperiodic on-site 
breather at the upper boundary of QC. The half-site breather exhibited a degeneration similar 
to that of the complex cyclic at the boundary to region H. The breather spontaneously changed 
back into an on-site breather and then decayed to rest. Multiple antisymmetric half-site breathers 
were also observed. A state with two half-sites 10-15 sites apart (center-to-center) were believed 
to be transferring energy by soliton shedding. A growth in amplitude in the center sites of one 
breather corresponded to the arrival of a soliton from the other breather and the expected 


decrease in the center-site amplitude of the second breather due to soliton ejection was observed. 


The upper and lower boundaries compare well with the theoretical boundaries (Fig. 
III.C.4). The upper boundary does not quite follow the hyperbola at the higher drive amplitudes, 
however the theory is approximate at these amplitudes. The lower boundary, which corresponds 
to the lattice decaying to rest, coincides almost exactly with the predicted threshold of lattice 
motion. Comparison of theoretical modulation envelopes and the numerical data were 
completed for extremely high, high and low response amplitude cases as well as a half-site case. 
The observed amplitudes were adjusted to that of the turning point as described in Section II.B. 
The normal 180° phase difference between lattice sites was eliminated for easier modulation 
envelope comparison. The extreme caSe was not as sharp as predicted and was = 70% of the 


theoretical amplitude (Fig. III.C.5). Similarly, the amplitude of the high case was ~90% of 
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theoretical but the envelope was sharper than expected (Fig. III.C.6). The low amplitude case 
exhibited the inverse characteristics of the extreme case (Fig. III.C.7). The observed envelope 
was sharper than the theory and =5% higher in amplitude. A comparison of the half-site case 


showed very good agreement with theory (Fig. III.C.8). 
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Figure III.C.6 Upper Cutoff Breather Theory Comparison 
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Figure III.C.8 Upper Cutoff Breather Theory Comparison 
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D. SOLITON SHEDDING 

The initial states for both upper and lower cutoff mode response investigations were 
evanescent states in a lattice of 100 sites. One end of the lattice was driven and the other end 
was given a free boundary condition. The large number of sites minimized reflection while 
maintaining the speed of the simulation. The parameters for the upper cutoff were y=.75, 
B=.03, while the lower cutoff were y=.50, 8=.03. Snapshots in time of the on-screen motion 
of the lattice were obtained to help describe the motion. In these "snapshot" figures, each dot 
is an individual lattice site and the site furthest to the left is the driven end site. 

Soliton shedding was observed in both the upper and lower cutoff, end driven lattices. The 
shedding was similar for both modes, with differences in the actual method of energy transfer 
from the driven site. The upper cutoff shedding was produced by the cyclic relaxation in 
amplitude of the driven end site. The lower cutoff shedding was also generated by end site 
relaxation but transferred energy through a "pivot" point before the soliton was created. 


An artificial potential well of the form: 


aC, nn cope 


iL. ‘pebee 


was used to replace f(x)= -x + x° in the lower cutoff program to permit amplitudes of the 
magnitude necessary to observe shedding in the softening mode. Otherwise, the oscillators would 
"go over the top." The upper threshold was marked by a distinct motion change and soliton 
Shedding (Fig. III.D.1). The shed solitons were very small and strongly evanescent from w= .95 
down to w=.6. The soliton shedding at w=.6 was distinct and the structure easily visible until 
20-30 sites from the driven end. The shedding results were similar for w= .5 with solitons visible 


until approximately site 35. A third harmonic component appeared at low drive amplitudes 
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(n= .144) and continued until the shedding threshold. The harmonic was identified by taking the 
FFT of a time series from the driven site (Fig. III.D.2). The visible result of the harmonic was 
an evanescent ripple wave (Fig. III.D.3). 

A fifth harmonic was identified for w=.3 (Fig. HI.D.4) and it also disappeared at the 
shedding threshold. The strength of the harmonic appeared to vary slowly with amplitude but 
over a fairly large range. The driven end site became quasiperiodic from »=.415 to 1.42 and 
from 2.14 to the shedding threshold for w=.2. Small wave packet ejection was observed from 
n=2.00 to 2.59. Upon reaching the shedding threshold, the end site relaxed and a large soliton 
was shed (Fig. III.D.5). The soliton would propagate 6-8 sites, then appear to shed a much 
smaller soliton. Both the original and shed soliton would then damp out quickly (Fig. III.D.6). 
A seventh harmonic appeared for w=.2 but was not visually identifiable. Odd harmonics up to 
the 11th were identified for w=.1 (Fig. I1I.D.7-8) and visually affected the lattice motion. The 
end site quasiperiodicity started at 7»=2.25 and the evanescent wave packets started at n=1.35. 
The shedding threshold was not reached by n=6.0. 

The dynamic threshold for the upper cutoff lattice exhibits a consistent trend even though 
the motion of the resultant state changes significantly between »=2.32 and 2.33 (Fig. III.D.9). 
The basic motion changes from soliton shedding to a non-shedding evanescent state. Shedding 
for w=2.025 occurs very slowly and the soliton propagates 30-40 lattice sites and then appears 
to stop (Fig. III.D.10). Apparently the nonlinearity lowers the frequency enough that linear wave 
propagation can occur inside the structure and it appears to "ring" as it decays (Fig. III.D.11). 
The shedding occurs slowly for w=2.05 and the soliton is large but jumbled (Fig. IIJ.D.12) and 
dies out quickly. The size of the shed solitons is cyclic for w=2.1-2.3 and the structure is very 
diffuse (Fig. III.D.13). Several small solitons are ejected followed by one large one and the 


relaxation of the end site is very distinct. Intermittent shedding occurs for w=2.31 and the first 
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8-9 sites move together (Fig. III.D.14). The transition point from shedding to nonshedding 
seems to be between w=2.32 and 2.33. The motion for w=2.32 is varied in the number of sites 
in tandem, the size of the shed solitons and the magnitude of the end site relaxation. The 
shedding is transitional for w=2.33 and the motion settles into a complex evanescent state. The 
first 6-9 sites move in tandem and the last site is the first of 2-3 evanescent sites. The motion 
for w=2.35-2.6 is similar but only 2-3 sites are in tandem. The vertical boundary between 
shedding and non-shedding states was confirmed by obtaining high drive level evanescent states 
and crossing the boundary into the shedding region. Specifically, the states w=2.35 7=5.05 and 
w=2.5 4=8.0 still exhibited the same lattice motion initiatated at the transition boundary. The 
States were then decremented in frequency to cross into the shedding region. The motion of both 


states slowly evolved to the previously described lattice motion for the final frequency. 


41 


+ — AZTIIBDIpotrusdisenb efaug x — UotrzI9? Sa yaxyIedanenm 
O — 3TQIstn sotruowuey 
6°6 3° 2°6 9°6 5° @ >’ G E°@ Z2°G | ae: 





SplLOYSeay, 2Or7zzeT Usntug—pug 


e349 


‘yFFoRND wanm07T 


Driven Lower Cutoff Drive Plane 





Figure III.D.1 


End- 


42 


Spectrum for clement @ LATTIC2X .C 
Alpha: 1.G0G600 Beta: 90.G30000 Gamma: © .500000 Eta: ©@.600000 Omega: G.5000080 
Max amp = 7069.307475 Max freq shown: 1.989437 Tiree interval: 0.125664 


Fundamental 





Figure III.D.2 Lower Cutoff FFT Spectrum w=.5 


Eta=.700 Omega=.500 





Figure III.D.3 Lower Cutoff Harmonic Ripple w=.5 
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Figure III.D.4 Lower Cutoff FFT Spectrum w=.3 
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Fta=3.7 Omega=.2 


Figure III.D.5 Time Snapshots w=.2 
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Figure III.D.6 Time Snapshots w=.2 
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Figure III.D.7 Lower Cutoff FFT Spectrum w=.1 


Eta=1.55 Onega=.1 





Figure III.D.8 Harmonic Ripples w=.1 
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Eta=.61 Omega=2.025 





Figure III.D.10 Time Snapshots w=2.025 
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Figure III.D.11 Time Snapshot w=2.025 
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Figure III.D.12 Time Snapshot w=2.05 
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Figure III.D.13 Time Snapshot w=2.2 
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Figure III.D.14 Time Snapshot w=2.31 


Eta=3.3 Omega=2.4 





Figure III.D.15 Time Snapshot w=2.4 
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IV. SUMMARY AND CONCLUSIONS 

It has been shown numerically that nonpropagating steady state breather solitons can occur 
in one-dimensional lattices of coupled nonlinear oscillators that are damped and parametrically 
driven. The simplicity of the model suggests that these states can exist in a variety of systems. 
Amplitude data agree well with a nonlinear Schrédinger theory at low amplitudes, but disagree 
substantially at higher amplitudes. The theory assumes that the amplitudes are weakly nonlinear 
and slowly varying in space. The continued existence of breathers at amplitudes where the theory 
breaks down is evidence of the robustness of these states. Dissipation and drive stabilize the 
breathers, although this is not yet understood. 

Breathers exist for a range of drive parameters (amplitude and frequency), and these 
regions have been mapped for both lower and upper cutoff breathers. Each region substantially 
disagreed with the theoretical prediction and, moreover, the natures of the two numerical regions 
were strongly dissimilar. The lower cutoff breathers are found to have a low-amplitude 
instability that is not predicted by the theory. The drive parameter region of the upper cutoff 
breathers contains a relatively large “finger” in which quasiperiodicity occurs, and a relatively 
large “island” in which the motion decays to rest. The corresponding instabilities are not yet 
understood, but it appears that the quasiperiodicity is a consequence of soliton shedding even 
though the shedding is not apparent except in a small region of the drive plane. This is based 
on the observation of strong quasiperiodic amplitude modulation during shedding and relatively 
weak modulation when shedding is not visible. Clearly visible soliton shedding is not expected 
to occur for a global parametric drive because a propagating breather has a spatially varying 


phase, whereas the drive is spatially constant. 
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For an end driven lattice, soliton shedding is dramatically visible. There is a drive 
amplitude threshold for the shedding to occur, and the value decreases as the linear frequency is 
approached. The phenomenon, although not yet understood, is expected to occur in any 
nonlinear wave system that possesses breathers. A possible application is the regeneration of 


fiber optic solitons which attenuate with distance. 
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APPENDIX 


PROGRAM COMMANDS 


I. 


2. 


a - Manual frequency increment, adjustable with A. 

A - Manual frequency increment adjustment, can be positive or negative. 

d - Coarse damping adjustment (decreasing), nominally .01, can be changed in 
variable declaration section of source code (the variable is Beta_increment). 

D - Coarse frequency adjustment (decreasing), nominally .001, can be changed in 
variable declaration section of source code (Omega _increment). 

Ctrl D - Coarse drive amplitude adjustment (decreasing), nominally .001, can be 
changed in variable declaration section of source code (Eta_increment). 

e - Fine damping adjustment (decreasing), 1/10 of coarse damping increment, can 
be changed in the mainQ section of the source code under the appropriate key hit 
section. 

E - Fine frequency adjustment (decreasing), 1/10 of coarse frequency increment, can 
be changed in the main() part of the source code under the appropriate key hit 
section. 

Ctri E - Fine drive amplitude adjustment (decreasing), 1/10 of coarse drive 
amplitude increment, can be changed in the mainQ section of the source code under 
the appropriate key hit section. 

Ctrl F - Time series record of one lattice site, prompts for a file name to store the 


series. This function works only in Text Mode (Ctrl T). The number and 


ay 


10. 


ee 


Ir 


Pee 


14. 


fa 


16. 


Ee 


20. 


periodicity of the samples are controlled in the eqn_motion( function section of the 
source code. Set for 80 samples, taken at the rate of one per time increment. 

G - Coupling adjustment (gamma), displays current value asks for new value. 
Ctrl G - Graphics mode, displays real time lattice motion. Also refreshes the screen 
without interfering with the lattice motion. Automatically invoked when the 
program is started or restarted. 

Ctrl H - Total lattice energy monitor, can only be used in Text Mode. 

i - Zoom in. The amplification of the screen is 2° for n button pushes. The 
amplification is vertical only. 

Ctrl I - Increases the strobing frequency of the phaseplot option, effectively 
increasing the sampling rate. 

Alt K - Coupling modulation, Options are random, gradient or sine wave. 

Ctrl K - Amplitude kick, specified lattice site by the desired amount. Amplitude and 
lattice site input from the keyboard. 

Ctrl L - Pins specified lattice site. Acts like a toggle switch, reselection unpins the 
Site. 

n - Random perturbation of the lattice. Perturbs all sites with a random coordinate 
change of < +3% of the maximum amplitude in the lattice. 

0 - Zoom out. The reduction of the screen is 2° for n button pushes. The reduction 
is vertical only. 

Ctrl O - Decreases the strobe frequency of the phaseplot option, effectively 


decreasing the sampling rate. 
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21. 


Dds 


23: 


24. 


25% 


26. 


ZT, 


28. 


Zo. 


30. 


31. 


p - Pause program. Freezes execution of the program. Alternate method of 
obtaining a time series file. 

P - Phaseplot mode. Selects up to five lattice site for phase space monitoring. 
Ctrl Q - Exits the program. 

r - User initiated state save. The lattice must be stable to within 1% (blue color) or 
the option will freeze the program until the 1% criterion is met. The option saves 
the current lattice modulation at the upper turning point of the lattice site being 
monitored. 

Ctrl R - Restarts the program without going back out to DOS. Goes to the same 
screen that was displayed when the program was first started up. 

s - Program freeze, captures lattice modulation when pushed. Any key will continue 
the program but the lattice will be at rest. 

S - Monitors the stability of the specified lattice site. Works like a toggle switch. 
The lattice will change colors when the monitored site is within 1% of its average 
amplitude for the last twenty upper turning points. 

Ctrl S - Turns off the drive. 

Ctrl T - Text Mode. Prints system parameters on screen. Total lattice energy can 
be monitored in this mode. 

u - Coarse damping adjustment (increasing), nominally .01, can be changed in 
variable declaration section of source code (variable is Beta_increment). 

Ctrl U - Coarse drive amplitude adjustment (increasing), nominally .001, can be 


changed in variable declaration section of source code (Eta_increment). 


a3 


32. 


33. 


34. 


35. 


36. 


ae 


38. 


39. 


40. 


41. 


42. 


U - Coarse Frequency adjustment (increasing), nominally .001, can be changed in 
variable declaration section of source code (Omega increment). 

v - Fine damping adjustment (increasing), 1/10 of coarse damping increment, can 
be changed in the main(Q) section of the source code under the appropriate key hit 
section. 

Ctrl V - Fine drive amplitude adjustment (increasing), 1/10 of coarse frequency 
increment, can be changed in the mainQ section of the source code under the 
appropriate key hit section. 

V - Fine drive amplitude adjustment (increasing), 1/10 of coarse increment, can be 
changed in the mainQ part of the su rce code under the appropriate key hit section. 
w - Displays peak lattice amplitude, works in Graphics Mode only. 

Ctrl W - Turns on waterfall type display. Trends are easily noticed on this type of 
display. 

Ctrl X - Selects lattice site to monitor for spectrum analysis. 

y - Dump spectrum generated by Ctrl X. 

z - Sets damping to zero. 

Z - Sets frequency to zero. 

Ctrl Z - Zeros the peak amplitude variable (variable automatically zeros when a 


parameter is changed). 
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B. PROGRAM CODE 

The QuickC program codes for the aforementioned numerical implementations (Sect III) 
are listed below. The full program listing for Lower Cutoff (CCLATL.C) is included and the 
equation of motion for Upper Cutoff (CCLATU.C), Lower Cutoff end-driven (ENDL.C) and 


Upper Cutoff end-driven (ENDU.C). 


1. Lower Cutoff, Global Drive 


/* PROGRAM LATTICE (VGA) 
VERSION 2.0 (QUICKC) 
#define LAST UPDATE 22 JUL 1991 BY CLEON WALDEN 
/* 
THIS PROGRAM SIMULATES A GENERAL LATTICE WITH EQUATIONS OF MOTION 
THAT CAN BE SUBSTITUTED IN WHERE INDICATED. VARIOUS INTERACTIVE 
FEATURES ARE PROVIDED, WHICH ARE EXPLAINED IN THE PROGRAMMER’S 
MANUAL FOR LATTICE. ENERGY MONITORING IS ADDED TO TEXT SCREEN. 
A WATERFALL DISPLAY IS ADDED TO MONITOR SOLITON MOTION DUE TO 
MEDIUM EFFECTS. 
*/ 


finclude "BIOS.H" 
#include "“graph.h" 
#include "stdlib.h" 
#include "CONIO.H" 
finclude "MATH.H" 
finclude "STDIO.H" 
#include "time.h" 
#include “direct.h" 
#include “process.h" 
#include “dos.h" 


#define getrandom(min,max) ((randQ % (int)((max)-(min)))+(min)+ 1) 
#define SQR(a) ((a)*(a)) 

#define CUB(a) ((a)*(a)*(a)) 

#define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr 

#define DOFOR(i,to) for(i=0;i <to;i+ +) 

#define PI 3.14159265359 

#define SCREEN CORRECTION FACTOR 1.05 

#define eta_increment 0.01 

#define beta_increment 0.01 
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#define omega _ increment 0.001 

#define STABILITY INCREMENT 0.01 
#define DISPLAY INCREMENT DINC 
#define PHASE INCREMENT PINC 
#define MIDDLE ELEMENT (no_pendulums/2) 
#define text_flag flags[0] 

#define graphics flag flags[1] 

#define file dump _flag flags[2] 

#define stability flag flags[3] 

#define peak flag flags[4] 

#define stable flag flags[5] 

#define phase flag flags[6] 

#define pause _flag flags[7] 

#define pause flag2 flags[8] 

#define stop flag flags[9] 

#define spectrum_flag flags[10] 

#define dump_spectrum_flag flags[{11] 
#define energy flag flags[12] 

#define waterfall flag flags[13] 

#define wait flag flags[14] 


/* The dynamical variables are entered here. */ 


double coordinate[150],momentum[150],old_coordinate[150],old_momentum[150], 
oldold_coordinate[150],oldold_momentum[150]; 


double acceleration, gamma[150],mean_gamma,beta,omega0[150],mean_omega0, 

eta,omega,alpha,model_time,time_int, 

max_amp,mode_amp,time series disp[8000],phase_elements[5], 

peak_record[20],pinned_elements[150], DINC,PINC, period, 

spectrum[4098] ,temp2[4096], energy, gradient,wat_inc, 

peak _amp,amp _ kick,freq_inc,delta; 

int no_pendulums,flags[15],counterl,phase_counter,phase_int,first_element, 

chosen_element,stability element,disp color,colors[5],counter2,step_size, 

spectrum _element,spectrum_counter,sample counter,energy_counter,g,gst[8000], 
waterfall counter,color_counter,number clicks,p flag,t_flag,top_flag; 


main( { 

FILE “fp; 

int c,i,j,k,l,m,n,xx,a; 

double modulation; 

char ans[5],buff[50],buff2[50]; 

/* void display textQ,stability_restartQ,user_initQ,display_graphicsQ, 
process_pause(),monitor_stability(,start_file dumpQ, 
stop_file_dumpQ,pin_elements(),kick_elements(Q),stability checkQ, 
init_phasplotQ,phasplotQ,compute_spectrum(),plot_spectrum() */ 
_clearscreen(_ GCLEARSCREEN); 
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user_init(); 

_setvideomode(_ MRES4COLOR); 

_setcolor(1); 

sprintf(buff,"Eta %.31f Omega %.3lf Gamma %.2I!f",eta,omega,mean_gamma); 
_registerfonts("TMSRB.FON"); 
_moveto(10, 10); 
_Setfont("t’tms rmn’bn5"); 
_outgtext(buff); 

xx =0; 

wat_inc=1; 

disp_color=1; 

stop_flag=0; 

counter2 =0; 

energy flag=0; 

energy _counter=0; 

color_counter=0; 

while(stop flag==0) { 

energy=0; 

if(pause flag2==1) { 

pause flag2=0; 

pause _flag=Q; 


process _pause(); 

} /* if(pause...) */ 

if(kbhitQ! =0) { /* Check for keystroke from user */ 
c=getchQ; 

if(c==17) { aa 0, > Quit program */ 

stop flag=1; 

else if(c= =112) { jean : Pause program/dump state */ 
pause flag=1; 

} 

else if(c==20) { oo : Text Mode */ 
_clearscreen(_ GCLEARSCREEN); 

text flag=1; 


graphics flag=0; 
waterfall flag=0; 
wait_flag=0; 
spectrum flag=0; 
display_textQ; 


} 

else if(Cc==7) { [* *G : Graphics Mode */ 
screen_graph(); 

phase flag=0; 

wait _flag=0; 

spectrum flag =0; 

waterfall flag=0; 

text_flag=0; 
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graphics flag=1; 
sprintf(buff, "Eta %.3lf Omega %.3lf Gamma %.2If" ,eta,omega,mean_gamma); 
_registerfonts("TMSRB.FON"); 


_moveto(10, 10); 

_setfont("t’tms rmn’bn5"); 

_outgtext(buff); 

} 

else if(C==8) { /* “H : Monitor energy in text mode */ 


if(energy_flag= =O) energy _flag=1; 
else energy flag=0; 


else if(C==16) { /**P : Kick in parameter variations */ 

Screen _textQ); 

srand((unsigned)time(NULL)); 

printf("\nEnter 1 if you wish to vary coupling: "); 

scanf(" %d" ,&j); 

if(==1) { 
printf("\nEnter 1 (random), 2 (gradient), or 3(sine) desired: "); 
scanf(" %d" ,&k); 


printf("\nEnter gradient (percentage over entire lattice): “); 
scanf(" %If" ,&gradient); 
DOFOR(i,no_pendulums) { 
gamma|[i] = gamma[i]+mean_gamma*i*gradient/(100*no_pendulums); 
} /* DOFORS/ 

+ /* ifG==2) */ 

else if(k= =3) { 

printf("\nEnter modulation amplitude (percentage): "); 

scanf(" %1f" ,&gradient); 

printf("\nEnter integer number of wavelengths to use: "); 
scanf(" %d" ,&l); 

DOFOR(i,no_ pendulums) { 
gamma|[i]= gamma[{i]+mean_gamma*gradient/100*(-cos(2*1*PI*1/ 
no_pendulums)); 


} 
else DOFOR(i,no_pendulums) { 
j=randQ; 
gamma[i]=mean_gamma*(.95+.1*j/RAND_ MAX); 
} 
} 


else { 
printf("\nVarying omegao...\n"); 
printf("\nEnter 1 (random), 2 (gradient), or 3(sine) desired: "); 
scanf(" %d" ,&k); 


printf("\nEnter gradient (percentage over entire lattice): "); 
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scanf(" %1f" ,&gradient); 
DOFOR(i,no_pendulums) { 
omega0{i] =omega0[i]+mean_omega0*i*gradient/(100*no_pendulums); 


} 


} 
else if(k= =3) { 
printf("\nEnter modulation amplitude (percentage): "); 
scanf("%1f" ,&gradient); 
printf("\nEnter integer number of wavelengths to use: "); 
scanf("%d",&1); 
DOFOR(i,no_pendulums) { 
omega0[i] =omega0[i] +mean_omega0*gradient/100*(-cos(2*1*PI*i/no_pendulums)); 


} 


else { 

DOFOR(i,no_pendulums) { 

j=randQ; 

printf("Random number is: %d RAND MAX tis %d\n",j,RAND MAX); 
omega0[i]=mean_omega0*(.95+.1*j//RAND MAX); 


} 
} 


screen_graph(); 


else if(C==23) { /* “W : Waterfall Display Mode */ 
if(waterfall flag==0O) { 
phase _ flag=0; 

graphics flag=0; 
text_flag=0; 
spectrum_flag=0; 
waterfall flag=1; 
init_waterfallQ; 

} 

else { 
if(wait_flag==1) { 
wait _flag=0; 
init_waterfallQ; 


else waterfall _flag=0; 


} 

else if(c==80) { [poe : Phase Plot Mode */ 
if(phase_flag= =O) init_phasplotQ; 

else { 

phase_flag=0; 


spectrum_flag=0; 
graphics flag=1; 
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waterfall flag=0; 
wait flag=0; 
screen _graphQ; 
po elecea) 


else if(c==6) { (oA : Time series to file */ 
if(file dump_flag= =O) start_file dumpQ; 

else { 

stop file dumpQ; 

file dump_flag=0; 


} 

else if(c==47) { [2G : Change gamma */ 
screen_textQ; 

printf("\nEnter new gamma, gamma is %If",mean_gamma); 
scanf(" %1f" ,&mean_ gamma); 

DOFOR(i,no_pendulums) { 

gamma[i]=mean_gamma; 


screen_graphQ); 
stability _restartQ); 


} 
else if(C==119) { /*w : Peak amplitude */ 
_moveto(210, 10); 
_setfont("t’tms rmn’bn5"); 
sprintf(buff2,"Amp %.4If",peak amp); 


_outgtext(buff2); 

else if(Cc==26) { (ees & : Reset peak_amp variable */ 
peak_amp=0.0; 

} 

else if(c==24) { (uae ; Calculate spectrum */ 


spectrum flag=1; 

screen_textQ; 

printf("\nEnter number of element to be analyzed: "); 
scanf("%d",&spectrum_ element); 

screen_graphQ); 


else if(C==21) { aes | 0 : Increase drive (coarse) */ 
Save _State(); 
eta=eta+eta_increment; 

stability _restartQ; 


else if(Cc==4) { a5 : Decrease drive (coarse) */ 
save State(); 
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eta=eta-eta_increment; 
stability _restartQ; 


else if(c==22) { oN, : Increase drive (fine) */ 
Save_state(); 
eta=eta+eta_increment/10; 

Stability restart(; 


else if(C==5) { (72a : Decrease drive (fine) */ 
Save state(); 
eta=eta-eta_increment/10; 

stability_restartQ; 


} 
else if(c==26) { [Rens : Turn off drive */ 
save state(); 
eta=0; 
stability restartQ; 
else if(c==85) { a : Increase frequency (coarse)*/ 


save State(); 
time_wait(); 
time_int=time_int*200/period; 
omega=omega+omega_ increment; 
period = 2*PI/omega; 
time_int=time_int*period/200; 
stability _restartQ); 


else if(C==68) { ib, : Decrease frequency (coarse)*/ 
save_state(); 

time_wait(); 

time_int=time_int*200/period; 

oOmega=omega-Omega increment; 

period =2*PI/omega; 

time_int=time_int*period/200; 

stability_restart(); 


else if(c==86) { aN. : Increase frequency (fine)*/ 
Save State(); 

time_wait(); 

time_int=time_int*200/period; 

omega=omega+omega_increment/10; 

period = 2*PI/omega; 

time_int=time_int*period/200; 

stability_restart(); 


else if(c==69) { f~ Ee : Decrease frequency (fine)*/ 
save_state(); 
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time_waitQ); 

time _int=time_int*200/period; 
omega=omega-omega_increment/10; 
period =2*PI/omega; 
time_int=time_int*period/200; 
Stability_restartQ; 


else if(c==90) { VEZ : Set frequency to zero */ 
save StateQ); 
omega=0; 
Stability restart(); 


else if(C==65) { /*°A : Set manual frequency increment*/ 
Screen _text(); 
printf("Manual frequency increment is %1lf\n",freq_inc); 
printf(" New increment:"); 
scanf(" %If" ,&freq_ inc); 
screen_graphQ); 


else if(Cc==97) { /*a : Manual Frequency increment*/ 
Save State(); 

time_wait(; 

time_int=time_int*200/period; 

omega=omega-+ freq _inc; 

period =2*PI/omega; 

time_int=time_int*period/200; 

Stability_restartQ; 


else if(C==117) { /*u : Increase damping (coarse)*/ 
Save StateQ; 

beta=beta+beta_increment; 

Stability restartQ; 


else if(c==100) { aad : Decrease damping (coarse)*/ 
Save State(); 

beta=beta-beta_increment; 

stability_restartQ; 


else if(c==118) { [*v : Increase damping (fine)*/ 
Save State(); 

beta=beta+beta_increment/10; 

Stability restartQ; 


else if(c==121) { aay, : Dump spectrum */ 
if(spectrum_flag==1) dump _ spectrum_flag=1; 


else if(C==101) { pre : Decrease damping (fine)*/ 
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Save_state(); 
beta=beta-beta_increment/10; 
stability restart(; 


} 

else if(C==122) { a7 : Turn off damping */ 
Save _State(); 

beta=0; 

Stability _restart(); 


else if(C==12) { [seca le : Pin elements */ 
pin_elementsQ; 
stability restartQ; 


j 

else if(C==11) { [* “K : Kick elements */ 
kick_elements(Q); 

Stability _restartQ; 


} 

else if(c==105) { [* j : Zoom in */ 
DINC=DINC/2; 

PINC =PINC*2; 

wat_inc=wat_inc*2; 


else if(Cc==111) { Ke ;: Zoom out */ 
DINC=DINC*2; 

PINC =PINC/2; 

wat_inc=wat_inc/2; 


else if(c==110) { 0 : Perturb system */ 
/* Maximum +/-3% of peak value */ 

linear_kickQ; 

screen_graphQ); 

stability restartQ); 


else if(C==15) { [* *O : Decrease strobe freq */ 
phase_int=phase_int-1; 


else if(C==9) { cai | : Increase strobe freq */ 
phase_int=phase_ int+1; 


else if(C==18) { /* *R : Restart program */ 

Screen text(); 

user_initQ; 

sprintf(buff,"Eta %.3lf Omega %.3lf Gamma %.2If",eta,omega,mean_gamma); 
_registerfonts("TMSRB.FON"); 

_moveto(10, 10); 

_setfont("t’tms rmn’bn5"); 

_outgtext(buff); 
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if(phase_flag==0) setvideomode( MRES4COLOR); 
else _setvideomode(_ VRES16COLOR); 


} 

else if(c==115) { ics : Stop everything as is */ 
DOFOR(i,no_pendulums) { 

coordinate{i] =0; 

momentum{i]=0; 

+ /* DOFOR */ 

stability_restartQ); 


} 

else if(c= =83) { i= S : Monitor stability */ 

monitor_stabilityQ; 

sprintf(buff, "Eta % .3lf Omega % .3lf Gamma %.2If",eta,omega,mean_ gamma); 
_registerfonts("TMSRB.FON"); 

_moveto(10, 10); 

_Setfont("t’tms rmn’bn5"); 


_outgtext(buff); 
} 


else if(C==114) { /* r : User initiated state save */ 
save stateuQ; 
} 


} — /* if(kbhit...) */ 

if(wait_flag==0) { 

eqn_motion(; 

} /* if(wait_flag */ 

} /* while(stop_ flag... */ 

_setvideomode@ DEFAULTMODE); 
printf("PROGRAM COMPLETE AT %lf",model_time); 
\ /* MAINO */ 


user init) { 

FILE *fq; 

int c,i,j,k,l,m; 

char answer[1],ans[{1],answ[1],filename[30]; 
printf("GENERALIZED LATTICE MODEL PROGRAM W/ VARIABLE PARAMETERS\n"); 
printf("Version 2.1X486 (MS) \n"); 

printf("Last updated 22 JUL 1991\n"); 

printf(" Variant notes: Default IC is AM\n"); 

printf(“Omega0 is set equal to one for all cases!\n"); 
printf("Rotating phase plane is used.\n"); 

printf("Real time FFT function is added...\n"); 

printf("Energy monitoring available via ~H in text mode\n"); 
printf("Set top flag to get output files for theory comparison. \n"); 
printf(" Manual Frequency increment is set to 1/20 omega inc.\n"); 
/* Sets flag to output additional information for theory comparison. 
Adds the information to the end of the normal output file*/ 
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printf("Set top_flag? "); 
scanf(" %s" ,ans); 
if((ans[0] = = 89) | | (ans[0]= =121)) { 
top_flag=1;} 
else top_flag=0; 
printf("\nDo you want to use a file for initial conditions (Y/N)? "); 
scanf(" %s" ,answer); 
sample_counter=0; 
if((answer[0] = =89)} | (answer[0]= =121)) { 
printf("\nOld file?"); 
scanf("%s" ,answ); 
printf(\nEnter name of file to be read: “); 
scanf("%s" ,filename); 

if((fq=fopen(filename, "r"))!=NULL) { 
mean_gamma=0; 
fscanf(fq," %d\n",&no_pendulums); 
fscanf(fq," %1f\n" ,&alpha); 
fscanf(fq," %If\n" ,&beta); 
fscanf(fq," %1f\n" ,&eta); 
fscanf(fq, %1f\n" ,&omega); 
DOFOR(i,no_pendulums) { 
fscanf(fq,"%lf %lf %if %lf\n", 
&omega0[i],&gamma|i],&coordinate[i],&momentum[i]); 
mean gamma=mean_gamma-+ gamma{i]; 
} 
mean gamma=mean_ gamma/no_pendulums; 
fclose(fq); 


else printf("Can’t open file requested."); 

} /* if((ans... */ 

else { 

printf("\nEnter number of pendulums to use: "); 
scanf("%d",&no_pendulums); 

printf("\nEnter mode amplitude: "); 

scanf(" %lf" ,&mode_amp); 

printf("\nEnter modulation amplitude: "); 

scanf(" %1f" ,&max_amp); 

printf("\nEnter nonlinear coefficient alpha (+/- 1 ONLY): "); 
scanf(" %1f" ,&alpha); 

DOFOR(k,no_pendulums) { 

if(alpha < 0) coordinate[k] =pow((-1),k)*(mode_amp 
+max_amp*sin(2*PI*k/no_pendulums)); 

else coordinate[k]=mode_amp+max_amp*sin(2*PI*k/no_pendulums); 
momentum|[k] =0; 

pinned elements[k]=0; 

} /* DOFOR */ 

printf("\nEnter coupling coefficient gamma: "); 
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scanf(" %1f" ,&mean_ gamma); 
DOFOR(i,no_pendulums) { 
gamma[i]=mean_gamma; 
omegaO[i]=1; 


printf("\nEnter drive amplitude eta: "); 
scanf(" %lf" ,&eta); 
printf("\nEnter drive frequency omega: "); 
scanf(" %1f" ,&omega); 
printf("\nEnter dissipation constant beta: “); 
scanf(" %1f" ,&beta); 
\ yemelse 3) 
printf("\nEnter number of steps per response cycle: "); 
scanf(" %d" ,&step_size); 
if ((step_size/2) != 0) g=.5; 
else g=0; 
time _int=2*PI/(step_size*omega); 
DOFOR(i, 15) flags[i] =0; 
DOFOR(i,4096) spectrum[i]=0; 
DOFOR(i, 150) { 
oldold_coordinate[i]=0.0; 
old_coordinate[i]=0.0; 
} 
screen_graph(); 

number_clicks=0; 

mean_omega0=1.0; 
peak_amp=0; 
spectrum_counter=0; 
/* Phase correction to start files in the proper phase*/ 
if((answ[{0] = = 89)] | (answ[0] = =121)) delta=asin((omega*beta)/eta)/2; 
else delta=(asin((omega*beta)/eta)/2)-PI; 
model_time=0.0; 
freq _inc=omega_increment/20; 
DINC = .02; /* CHANGE THIS TO CHANGE SCALE OF DISPLAY */ 
PINC=2; /* CHANGE THIS TO CHANGE SCALE OF PHASE PLOT */ 
} /* USER INIT */ 


display graphicsQ { 
/* Displays the lattice sites. The screen updates sequentially through the 
points after each round of motion calculations are complete. */ 
int c,i,j,k,l,m,n; 
if(no_pendulums < =40) { 
first_element=0; 
DOFOR(k,no_pendulums) { 
n=0; 
if((stability_element= =k)&&(stability flag==1)) n=1; 
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if((pinned_elements[k]= = 1)&&(disp_color==1)) n=2; 
if((pinned_elements[k]= = 1)&&(disp_color= =2)) n=-1; 
1=old_coordinate{k]/DISPLAY_INCREMENT; 
_setcolor(0); 
_setpixel((160-S*MIDDLE ELEMENT + 5*k),(100+])); 
1=coordinate{k]/DISPLAY_INCREMENT; 
_setcolor(disp_color+n); 
_setpixel((160-S*MIDDLE ELEMENT + 5*k),(100+1)); 
} /* DOFOR */ 
elk */ 

else { 
|=(no_pendulums-40)/2; 
first_element=]; 
DOFOR(k,40) { 
m=old_coordinate{l+k]/DISPLAY INCREMENT; 
n=Q; 
if((stability element= =(1+k))&&(stability flag==1)) n=1; 
if((pinned_elements[{k]= = 1)&&(disp_color==1)) n=2; 
if((pinned_elements[{k] = = 1)&&(disp_color==2)) n=-1; 
_setcolor(0); 
_setpixel((60 + 5*k),(100+m)); 
m=coordinate{l+k]/DISPLAY_ INCREMENT; 
_setcolor(disp_color+n); 
_setpixel((60 + 5*k),(100+m)); 

} /* DOFOR */ 
eee LSE */ 
} /* DISPLAY GRAPHICS */ 


display_textQ { 

/* The coordinate and other information below is a snapshot of what the 
lattice is like at the particular instant in time when the key was pushed.*/ 
char message[80]; 

int c,j,k,1,m; 

_setvideomodel DEFAULTMODE); 

printf("Time is : %lf",model_time); 

printf(" System parameters are: \n"); 


printf(" Gamma %\f" ,mean_gamma); 
printf(" Eta %\f" eta); 

printf(" Omega % \f\n" ,omega); 
printf(" Beta %\f" beta); 

printf(" Alpha %\f\n" ,alpha); 


printf("\nThere are %d elements in the system" ,no_pendulums); 
printf("\n\nPress any key to continue...”); 

c=getchar(); 

while(kbhitQ = =0) ; 


printf("Element Position Velocity Element Position Velocity\n"); 


DOFOR(j,20) { 
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printf(" %d %olf %olf %od %olf %\f\n" , 
j,coordinate[first_element+j}],momentum[first_element+j}, 

(j+ 20),coordinate[first_element+j +20], 

momentum|[first_element +j+20]); 

} /* DOFOR */ 

printf("\n\nPress ~G for graphics, ~Q to quit."); 

c= getchar(); 

while(kbhit0 = =0); 

while(kbhitQ = =0); 

\ /* display_text */ 


process pause() { 
/* Old method of saving files. Not sure if it has other uses.*/ 
int c,i,j,k,l,mode; 
FILE *fr; 
char filename[30}], message[80]; 
_clearscreen@ GCLEARSCREEN); 
_setvideomodel DEFAULTMODE); 
printf("Do you wish to save this state? "); 
if((c=getcheQ) = =121) { 
printf("\nEnter name of file to be written: "); 
scanf("%s" filename); 
if((fr=fopen(filename,"w"))!=NULL) { 
fprintf(fr," %d\n" ,no_pendulums); 
fprintf(fr, " %1f\n" ,alpha); 
fprintf(fr," %1f\n" ,beta); 
fprintf(fr, " %1f\n" eta); 
fprintf(fr, " %1f\n" ,omega); 
DOFOR(i,no_pendulums) 
fprintf(fr,"%lf %1f %lf %lf\n", 
omega0[i],gammafi],old_coordinate[i], momentum[i]); 
fclose(fr); 
EZ) 3) 
else printf("Failed to open %s\n",filename); 
ps1) % 
time_int=time_int*200/period; 
printf("\nEnter new time multiple (old multiple is %1f): ",time_int); 
scanf(" %1f" ,&time_int); 
time _int=time_int*period/200; 
if(phase flag==0) setvideomode@ MRES4COLOR); 
else _setvideomode( VRESI6COLOR); 
\ /* process pause */ 


start file dumpQ { 

/* Sets up a file dump of the coordinates of the selected element. Set up 
to record every time through for 30 times */ 

_clearscreen@ GCLEARSCREEN); 
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_setvideomode@ DEFAULTMODE); 
printf("Enter number of element to be monitored: "); 
scanf("%d" ,&chosen_ element); 
if(chosen_element > (no_pendulums-1)) 
printf("Out of range. No file dump"); 
else file dump_flag=1; 
_clearscreen(_ GCLEARSCREEN); 
if(phase_flag==0) 

_Setvideomode(_ MRES4COLOR); 

else 

_Setvideomode(_ VRES16COLOR); 
counter 1 =0; 

} /* start file dump */ 


stop file dumpQ { 

/* Dumps the collected data to a file. */ 

char filename[30]; 

int c; 

FILE *fs; 

_clearscreen(_ GCLEARSCREEN); 
_Setvideomodel DEFAULTMODE); 
printf("\nEnter name of file to be written: "); 
scanf(" %s" , filename); 

if((fs = fopen(filename,"w"))!=NULL) { 
DOFOR(c,8000) fprintf(fs," %If %df\n" time series disp[c],gst[c]); 
_clearscreen( GCLEARSCREEN); 

counter 1 =Q; 

if(phase_flag= =0Q) 

_Setvideomode(_ MRES4COLOR); 

else _setvideomode(_ VRES16COLOR); 

} /* if(Q) */ 

} /* stop_file dump */ 


monitor_stabilityQ { 

/* Selects the lattice site to be monitored */ 
if(stability flag==0) { 

stability flag=1; 

disp_color=2; 

_setvideomode( DEFAULTMODE); 

printf("Enter number of element to be monitored: "); 
scanf(" %d" ,&stability_element); 

Stability restart (); 


else { 
Stability flag=0; 
disp_color=1; 


} 
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if(phase_flag= =O) 
_setvideomode(_ MRES4COLOR); 
else 

_setvideomode(_ VRES16COLOR); 
j 


stability checkQ { 

/* The criterion is set for 1%. Appears to check to see if 20 consecutive 
peaks are within the criterion of each other. */ 

int 1,k; 

k=0; 

DOFOR(i,19) { 

stable flag=1; 

if((peak_record[i+ 1] >(peak_record[i]*(1 +STABILITY_INCREMENT)))! 
(peak_record[i+ 1] < (peak_record[i]*(1 - STABILITY INCREMENT )))) { 
stable flag =0; 

break; 

loint */ 

\ /* DOFOR */ 

if(stable flag==1) { 

disp _color=1; 

j 

} 


pin_elements() { 

/* Pins an element at rest */ 

int ans,i,check; 

_setvideomodel DEFAULTMODE); 

check=0; 

printf("\nEnter number of element to be pinned: "); 
scanf(" %d" ,&ans); 

if((ans <0)|(ans > (no_pendulums-1))) { 

check= 1; 


} 
if(check= =0) { 
if(pinned elements[ans]= =O) { 
pinned _elements[ans]= 1; 
coordinate[ans]=0; 
momentum|[ans]=0; 


else pinned_elements[ans]=0; 

ea) 

if(phase_flag= =O) 
_setvideomode( MRES4COLOR); 
else 

_setvideomode(_ VRES16COLOR); 


} 


70 


kick elementsQ) { 

/* Equivalent to banging an element with a hammer. No control over when 
the hit occurs in the cycle. */ 

int ans,check; 

double amount; 

_setvideomode@ DEFAULTMODE); 

check =Q; 

printf("\nEnter number of element to kick: “"); 
scanf(" %d" ,&ans); 

if((ans <Q)| (ans >(no_pendulums-1))) { 
check=1; 


} 

if(check= =0) { 

printf("\n Enter amount to kick: °); 
scanf(" %1f" ,&amount); 
coordinate[ans}=coordinate[ans] + amount; 
ha if.*/ 

if(phase_flag= =0) 

_setvideomode@ MRES4COLOR); 

else 

_Setvideomode(_VRES16COLOR); 


} 


stability restartQ { 
/* Changes the color to show restart and zeros the stability check matrix */ 
int 1; 
char buff1[{50}; 
disp _color=2; 
stable flag =0; 
DOFOR(i,20) peak_record{i]=0; 
peak_record[15]=10; 
peak_amp=0; 
sprintf(buffl,"Eta %.3lf Omega %.3lf Gamma %.2I1f",eta,omega,mean_gamma); 
_moveto(10,10); 
_Setfont("t’tms rmn’bn5"); 
_outgtext(buff1); 


save state() { 
/* Automatic state saver. It’s supposed to keep the program from crashing 
when a lattice flys off numerically. */ 
int 1; 
char filename[}]= "savedsta"; 
BILE *fl; 
/* change the number below to change strokes between autosaves */ 
if(number_clicks == 100) { 
if(chdir("E:\\MODEL\\CLEON") != 0) { 
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screen _text(); 

printf("failed to change to E:\\MODEL\\CLEON.\n"); 
stability restartQ; } 

p_flag=0; 

while(p_flag==0) { 

eqn_motionQ; 


} 
if((fl=fopen(filename,"w")) != NULL) { 
fprintf(fl," %d\n",no_pendulums); 
fprintf(fl, “ %1f\n" ,alpha); 
fprintf(fl," %1f\n" ,beta); 
fprintf(fl," %1f\n" ,eta); 
fprintf(fl," %1f\n" ,omega); 
DOFOR(i,no_pendulums) 
fprintf(fl,"%lf %If %lf M%lf\n", 
omega0[i],gammaf[i],old_coordinate[1],old_momentum[I]); 
fclose(fl); 
number clicks=0; } 
else { screen text(; 
printf("Failed to open %s\n",filename); 
stability restart; } } 
if(chdir("E:\\MODEL") != 0) { 
screen_text(); 
printf("failed to change back to E:\\MODEL.\n"); 
stability _restartQ; } 
else + +number clicks; 
screen graphQ; } 
/* save_state */ 


save StateuQ { 
/* User initiated state save. Normal way to save lattice data files. Saves 
data for the lattice at a turning point (within one time step). */ 
int 1; 
char filename[30}; 
PILE *fp: 
screen_text(); 
printf("\nEnter name of file to save state in: "); 
scanf(" %s" filename); 
p_flag =0; 
while(p_flag==0) { 
eqn_motionQ); 


} 
printf("\n%1lf %1f" ,model_time,delta); 
if((fp=fopen(filename,"w")) != NULL) { 
fprintf(fp,"%d\n",no_pendulums); 
fprintf(fp," %1f\n" ,alpha); 
fprintf(fp, " %1f\n" ,beta); 
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fprintf(fp," %1f\n" ,eta); 
fprintf(fp, " %1f\n" ,omega); 
DOFOR(i,no_pendulums) 
fprintf(fp,"%lf lf %lf Mlf\n", 

omega0[{i],gammaf[i],old_coordinate{i],old_momentun{Ii]); 

if(top flag= =1){ 

fprintf(fp,"%If %1lf %Ilf %lf\n",oldold_coordinate[stability element], 
old_coordinate[stability element],coordinate{stability element],time_int); 

fprintf(fp,"%lf %1f %I|f",oldold_momentum{stability element], 
old _momentun({stability element],momentum{stability_element]);} 

fclose(fp); } 

else printf("Failed to open %s\n",filename); 
pause flag=0; 
printf("\nEnter number of steps per response cycle: “); 
scant(" %d" ,&step_size); 
if ((step_size%2) != 0) g=.5; 
else g=0; 
time_int=2*PI/(step_size*omega); 
screen graphQ); } 

/* save stateu */ 


retrieve stateQ) { 


/* Retrieves the auto saved file when the lattice crashes. */ 
int 1,C; 
char filename[]= "savedsta"; 
HILE *fg; 
screen_text(); 
printf("\nLATTICE CRASH!!!!\n"); 
printf("Do you want to retrieve the latest stored state(Y/N).\n"); 
printf(" Eta %lf Omega %lf\n",eta, omega); 
if((c=getcheQ))==121) { 
if(chdir("E:\\MODEL\\CLEON") != 0) { 
printf("Failed to change directories.\n"); 


} 


else { printf("\nEnter name of file to retrieve.: "); 
scanf("%s" filename); } 

if((fg=fopen(filename,"r"))!=NULL) { 
mean_gamma=0; 

fscant(fg," %d\n",&no_pendulums); 

fscanf(fg," %1f\n" ,&alpha); 

fscanf(fg," %lf\n" ,&beta); 

fscanf(fg," %1f\n" ,&eta); 

fscanf(fg," %1lf\n" ,&omega); 

DOFOR(i,no_pendulums) { 

fscanf(fg,"%lf lf %Iif %lf\n", 
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&omega0{[i],&gamma|{i], &coordinate[i] ,&momentum{i]); 
mean_gamma=mean_gamma+ gamma{i]; 


fclose(fg); 

counter2 =0; 

delta = atan((omega*beta)/sqrt(SQR(eta)-(SQR(omega)*SQR(beta))))/2; 
model_time=0.0; 


else printf("Can’t open file requested."); 
time_int=time_int*200/period; 

printf("\nEnter new time multiple (old multiple is %If): “,time_int); 
scanf(" %1f" ,&time_int); 

time_int=time_int*period/200; 

if(chdir("E:\\MODEL") != 0) { 

screen_textQ; 

printf("failed to change back to E:\\MODEL.\n"); } 

screen_graphQ); 


/* retrieve state */ 


screen _textQ) { 
_clearscreen( GCLEARSCREEN); 
text_flag=1; 
graphics flag=0; 
spectrum_flag=0; 
_setvideomodel DEFAULTMODE); } 
/* screen text */ 


screen graphQ) { 
_clearscreen(_ GCLEARSCREEN); 
_setvideomode(_ MRES4COLOR); 
_setcolor(1); 
phase _flag=0; 
spectrum _flag=0; 
text_flag=0; 
graphics flag=1; } 
/* screen graph */ 


linear_kickQ { 
/* Lattice perturbation up to a maximum of 3%. Sort of random as to sign of 
kick direction for each site. Uses the % of the maximum amplitude site. */ 
int a,j,1; 
double k,b; 
struct dostime_t time2; 
_dos_gettime(&time2); 
srand(time2.hsecond); 
DOFOR(i,no_pendulums) { 
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b=randQ; 

/* change 33.33 to change %, now 3% */ 
k=b/(33.33 * 32767); 

amp _kick=k*peak_amp; 

if(o > 16384) { 

coordinate[i] = coordinate[1] +amp_kick; 


else coordinate[i} = coordinate[i}] - amp_kick; 


peak_amp=0.0; 
} /* linear_kick */ 


time_wait() { 
/* part of a routine to get the lattice at the turning point. */ 
t_flag=0; 
while(t_flag= =0){ 
eqn_motionQ; 


/* time wait */ 


init phasplotQ { 

/* Sets up the phase plot to look at the phase plane of sites. */ 
int c,i,j,k,1; 

char ans[5],message[40],txt[3]; 

float dummy; 


spectrum_flag=0; 

_setvideomodel DEFAULTMODE); 
DOFOR(i,5) phase_elements[i]=0; 
printf("\nDo you want Poincare sections? "); 
if((c=getch())= = 121) phase _flag=2; 

else phase flag=1; 


k=0; 

while((it+ +!=999)&&(k <5)) { 
scanf(" %d" ,&1); 
phase_elements[k+ +]=i+1; 

+ /* while */ 

DOFOR(i,5) { 
if(phase_elements[i] = = 1000) { 
phase_elements[i}=0; 

} 


} 
_Setvideomode(_ VRES16COLOR); 


_clearscreen(_ GCLEARSCREEN); 
graphics _flag=0; 
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printf("\nWhich elements do you wish to monitor (999 to finish): 


_settextcolor(7); 
colors(0] =2; 
colors[1]=3; 
colors(2}=9; 
colors[3]= 14; 
colors[4]= 13; 
DOFOR(i,120) { 
_setpixel(318,4*1); 


} 

DOFOR(i, 160) { 
_setpixel(4*i,238); 
} 


DOFOR(i,6) { 
_setpixel(319,10*SCREEN CORRECTION FACTOR+1 
+(80/SCREEN CORRECTION FACTOR)*(i+ 1); 


} 
DOFOR(i,8) { 
_setpixel(80*(i+ 1),239); 


j 
_moveto(480,450); 
DOFOR(i,5) { 
if((phase_elements[i])!=0) { 
c=phase_elements[i]-1; 

/[* itoa(c,txt, 10); 
_outtext(txt); */ 
printf("%d ",c); 
_Setcolor(colors[i]}); 
DOFOR(,4) { 
DOFOR(k,8) { 
_setpixel(550+ 10*i+k,10+j); 


wey ye Wey 


} /* DOFOR */ 

printf(" LAT2X486.C"); 

printf("\nGamma: %lf Eta: %lf Beta: %lf",mean_gamma,eta,beta); 
printf("\nOmega: %lf Alpha: %lf Time interval: lf", 
omega,alpha,time_ int); 

dummy = 1/(time_int*omega); 

phase _int=dummy/1; 

phase_counter=0; 


} 


phasplotQ { 

/* Displays the phase for chosen elements. */ 

int 1,j,k,1; 

double speed, position,fast_coordinate,fast_ momentum; 
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if(phase flag==1) { 

DOFOR(i,5) { 

if((j=phase_elements[i]-1)!=-1) { 

speed =PHASE INCREMENT*momentum|j]; 
_setcolor(colors[i]); 

position=PHASE INCREMENT *coordinate[j]; 
fast_ coordinate = (position*cos(omega*model_time) 
-speed*sin(omega*model_time)/omega) 
/(3*SCREEN CORRECTION FACTOR); 

fast_ momentum = (position*sin(omega*model_ time) 
+speed*cos(omega*model_time)/omega)/4; 
_setpixel((320+320*fast_coordinate),(240+240*fast_momentum)); 
} /* if(O) */ 

} = /* DOFOR */ 

of it */ 

else if(phase_flag==2) { 
if(phase_counter==phase int) { 
phase_counter=0; 

DOFOR(;,5) { 

if((j=phase_ elements[i]-1)!=-1) { 

speed = PHASE INCREMENT*momentum[}j]; 
_setcolor(colors[i]); 
position=PHASE INCREMENT*coordinate[j]; 
fast_coordinate= (position*cos(omega*model_time) 
-speed*sin(omega*model_time)/omega); 

fast_ momentum = (position*sin(omega*model_ time) 
+speed*cos(omega*model_time)/omega)/4; 
_setpixel((320+320*fast_coordinate),(240+240*fast_momentum)); 
} /* if() */ 

} /* DOFOR */ 

baie it */ 

} /* else ifQ */ 

} /* phasplot */ 


/* FFT Routines for lattice programs */ 


compute spectrum({ /* Uses algorithm from p. 411 of Numerical Recipes in C */ 
int n,mmax,m,j,istep,i,nn,temp; 

double wtemp,wr,wpr, wpi, wi,theta,tempr,tempi; 
nn=2048; 

n=nn< <1; 

io, 

DOFOR(i,2048) { 

nn=(1&1024)/1024; 

nn=nn+ (1&512)/256; 

nn=nn+ (i&256)/64; 

nn=nn-+ (1&128)/16; 


| 


nn=nn + (i&64)/4; 

nn=nn + (i&32); 

nn=nn+ (i&16)*4; 

nn=nn+ (i&8)* 16; 

nn=nn+ (i&4)*64; 

nn=nn+ (i&2)*256; 

nn=nn+ (i&1)*1024; 
temp2[2*nn]=spectrum[2*i]; 
temp2[2*nn+ 1] =spectrum[2*i + 1); 


} 

DOFOR(i,4096) spectrum[1i] =temp2[1}; 
mmax =2; 

while(n > mmax) { 

istep = 2*mmax; 

theta=2*PI/mmax; 

wtemp =sin(0.5*theta); 

wpr =-2.0*wtemp*wtemp; 

wpi=sin(theta); 

wr= 1.0; 

wi=0.0; 

for(m=0;m < (mmax-1);m+ =2) { 
for(i=m;i< =n;i+ =istep) { 

j=1+mmax; 
tempr=wr*spectrum{[j]-wi*spectrum|[j + 1]; 
tempi=wr*spectrum[j + 1] + wi*spectrum{[j]; 
spectrum[j] =spectrum[i]-tempr; 

spectrum{j + 1] =spectrum[i+ 1]-temp1; 
spectrum[1i] + =tempr; 

spectrum[i+ 1]+ =tempi; 


wr =(wtemp = wr)*wpr-wi*wpi+ wr; 
W1=wi*wpr + wtemp*wpi+ wl; 

° 

mmax = istep; 

j 

} 


plot_spectrumQ) { 

int 1,}; 

double freq,magnitude, mag ,max; 

int C; 

FILE *fp; 

if(dump_spectrum_flag==1) { 
printf("\nDumping spectrum now"); 

fp =fopen("spectrum.out" ,"w"); 
for(i=1;1<2048;1+ +) { 

if(time_int! =0) freq=i/(2048*time_int); 
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else printf("Time_int was zero!"); 

magnitude =sqrt(spectrum[(2*i)]*spectrum[(2*i)] 
+ spectrum[(2*i+ 1)]*spectrum[(2*i+ 1)]); 
fprintf(fp,"%1f %1f \n",freq,magnitude); 


} 
fclose(fp); 
dump _ spectrum _flag=0; 


} 

_clearscreenl GCLEARSCREEN); 

max =Q; 

for(i= 1;i1< 2049;i+ +) { /* this loop calculates the max spectral component*/ 
magnitude=sqrt(spectrum[(2*(i))]*spectrum[(2*(i))] 

+spectrum[(2*(i) + 1)]*spectrum[(2*(i) + 1)])); 

if(magnitude > max) max = magnitude; 


_setcolor(3); 

DOFOR(i,640) _setpixel(i,460); 

DOFOR(i,11) _setpixel(51*i,461), setpixel(51*i,462); 
printf("Spectrum for element %d LATTIC2X.C",spectrum element); 
printf("\nAlpha: %|lf Beta: %lf Gamma: %lf Eta: %1lf Omega: %If", 
alpha,beta,mean_gamma,eta,omega); 

printf(\nMax amp = %lf = ",max); 

freq=512/(2048*time_int); 

printf("Max freq shown: %lf Time interval: %If",freq,time_int); 
_setcolor(1); 

DOFOR(,512) { 

alll 

magnitude= sqrt(spectrum[2*i]*spectrum[2*1] 

+ spectrum[2*i + 1]*spectrum[2*i+ 1]); 

mag = magnitude*400/max; 

while(j + + < mag) { 

_setpixel(i,460-}); 

} /* while */ 

} /* DOFOR */ 

} 


init_waterfallQ { 

/* Sets up waterfall type display to see trends in latice motion. Good for 
observing kink or breather drift. */ 

int c,1,j,k,1,m; 

char ans[5]; 

_SetvideomodeC VRES16COLOR); 

_setcolor(1); 

DOFOR(i,600) _setpixel(20+i,20); 

DOFOR(i,450) _ setpixel(20,i+ 20); 

DOFOR(,61) { 

if((i= = 10)| G@= =20)] (i= =32)| (= =42)| G@= =52)) _ setcolor(2); 
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DOFOR(j,113) { 
setpixel(20 + 10*1,20+ 4*)); 


} 


_setcolor(1); 

} 

_setcolor(2); 

DOFOR(i,11) { 

DOFOR(, 160) { 
_setpixel(20 + 4*j, 16+ 40*(i+ 1)); 
} 


waterfall counter =0; 
color_counter=2; 


} 


disp_waterfallQ { 

int c,i,j,k,1; 
if(waterfall_counter+ + < 110) { 
if(color_counter + + = =14) { 
color_counter=2; 
_setcolor(color_counter); 

} 

else { 

_setcolor(color_counter); 


} 
DOFOR(i,no_pendulums) { 
if(waterfall counter <55) { 
_setpixel(20 + 2*1,24+ 8*waterfall_counter-7*wat_inc*(coordinate[i]})); 
} 
else { 
setpixel(340 + 2*1,24 + 8*(waterfall_counter-55)-7*wat_inc*(coordinate[i]})); 


aw ew =~ | 


else { 

wait flag=1; 
} 

} 


eqn_motionQ { 

/* The following loop calculates a round of velocities for one 
iteration of the model clock. This is where the equations 

of motion need to be located, if one wishes to change the 
model’s physics! */ 

int 1,j,k; 

DOFOR(i,no_pendulums) 

if(coordinate[i] > 100) { 


80 


retrieve state(); 


DOFOR(i,no_pendulums) { 
if(i= =0) { 
acceleration= gamma|i]*(coordinate[i+ 1]+coordinate[no_pendulums-1] 
-2*coordinate[i])-beta*momentum[i] 
-(SQR(omega0[i]) +2*eta*cos(2*omega*model_time-delta))*coordinate[i] 
+ alpha*CUB(coordinate[i]); 


} 

else if(i= =(no_pendulums-1)) { 

acceleration = gamma|i]*(coordinate[i-1]+ coordinate[0] 
-2*coordinate[i])-beta*momentum[i] 

-(SQR(omega0[i]) + 2*eta*cos(2*omega*model_time-delta))*coordinate[1] 
+ alpha*CUB(coordinate[1]); 


else { 

acceleration= gamma|{i]*(coordinate[1+ 1]+coordinate[i-1] 
-2*coordinate[i])-beta*momentum[1] 

-(SQR(omega0[i]) + 2*eta*cos(2*omega*model_time-delta))*coordinate[i] 
+ alpha*CUB(coordinate[1]); 


oldold_momentun{i]=old_momentun{i]; 
old_momentum[i]=momentun({i]; 

momentum{i] =momentum[i] + acceleration*time_int; 
\ /* DOFOR */ 

DOFOR(i,no_pendulums) { 
oldold_coordinate[{i]=old_coordinate[i]; 
old_coordinate[i] = coordinate[i]; 
if(pinned_elements[i]= =0) { 
coordinate[i] =coordinate[i] +momentum[i]*time_ int; 


} 


/* Records data for file dump */ 

if((i==chosen_ element)&&(file dump flag==1)) { 
/* if(counter2 + + = =30) { */ 

time series disp[counter1 + + ]=coordinate[i]; 

gst{counter! + +]=g; 

/* counter2=0; 

D/ 

if(counter 1 = =8000) stop file dumpQ; 

peti... */ 


/* Stores energy information */ 

if(energy flag==1) { 

energy =energy + 0.5*(SQR(momentun(i]) 

+ gamma|[i]*SQR((coordinate[i+ 1]-coordinate[i])) 
+ SQR(coordinate[i}))-alpha/4*pow(coordinate[i],4); 
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} /* if(energy...) */ 


/* Stores FFT information */ 

if((spectrum_flag= = 1)&&(spectrum_element= =i) 
&&((sample_counter+ +)==1)) { 
if(spectrum_counter < 2048) { 
spectrum[2*spectrum_counter] =coordinate[1]; 
spectrum[2*spectrum_ counter + 1]=0; 
spectrum_counter+ +; 

sample _counter=Q; 


else { 

compute _spectrum(); 
_clearscreen( GCLEARSCREEN); 
_setvideomode(_VRES16COLOR); 
plot_spectrumQ); 

phase_flag=0; 

text_flag=0O; 

graphics flag=0; 

DOFORG,4096) spectrum{[j] =0; 
spectrum counter =0; 
sample_counter =0; 

jo) aeise 

tal* it (Spec... */ 


/* Records the peaks for the stability checker */ 

if(((i= =stability element)&&(stability flag==1)&&(stable_flag==0))| 
((i= =(no_pendulums-1))&&(waterfall flag==1))) { 
if((coordinate[i] > 0)&&(coordinate[i] < old_coordinate[i]) 
&&(peak flag==0)) { 

peak flag=1; 

if(pause_flag==1) pause_flag2=1; 

DOFOR(k, 19) peak_record[k]=peak_record[k + 1]; 
peak_record[19]=coordinate[i]; 

if(stability_flag==1) stability _checkQ; 
if(waterfall_flag==1) disp_waterfallQ; 


else if(coordinate[i] <0) peak_flag=0; 


} 


/* Records peak amplitude for display and linear kick */ 
if((coordinate[i] > 0)&&(coordinate[i] < old_coordinate[i])) { 
if(old_coordinate[i] > peak_amp){ 

peak_amp = old_coordinate[i]; 


} 
} 
if((i= =stability element)&&(stability flag= = 1)&&(stable_flag==1)){ 
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if(((coordinate[i]-old_coordinatef[i]) > = 0) && 
((old_coordinate{i]-oldold_coordinate[i}) < = 0) && (coordinate{i] < 0)){ 
p_flag=1; 


} 
+ /* DOFOR */ 


/* Increments time and resets it at 2Pi. Causes problems in eqns of 
motion if not reset */ 

model_time=model_time+time_int; 

if((g+1) == (step _size/2)){ 

model_time=0.0; 


t flag=1; 
eae 0; 

} 

else o++; 


/* Displays energy information gathered above. */ 

if(graphics flag==1) display_graphicsQ; 
if((text_flag==1)&&(energy flag==1)) { 

printf("\n Total Energy at time %lf is %1f",model_time,energy); 
energy counter+ +; 

if(energy_counter==20) { 

printf("\n\nStrike any key to continue..."); 

while(kbhitQ = =0Q); 

energy counter=0; 


} 

if(phase_ flag!=0) { 

phasplotQ; 

phase _counter++; /* Used for Poincare sections */ 
} 

} 


2. Upper Cutoff, Global Drive, Equation of Motion 


eqn_motionQ { 
/* The following loop calculates a round of velocities for one 
iteration of the model clock. This is where the equations 
of motion need to be located, if one wishes to change the 
model’s physics! */ 
int i,j,k; 
DOFORG,no_pendulums) 
if(coordinate[i] > 100) { 
retrieve stateQ; 
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DOFOR(i,no_pendulums) { 
if(i= =O) { 
acceleration = gamma|[i]*(coordinate[i+ 1]+coordinate[no_pendulums-1] 
-2*coordinate[i])-beta*momentum[i] 
-(SQR(omega0[1}) + 2*eta*cos(2*omega*model_time-delta))*coordinate[i] 
+ alpha*CUB(coordinate[1]); 


else if(i= =(no_pendulums-1)) { 

acceleration = gamma|i]*(coordinate[i-1]+ coordinate[0] 
-2*coordinate[i])-beta*momentum|i] 

-(SQR(omega0[i]) +2 *eta*cos(2*omega*model_time-delta))*coordinate[i] 
+ alpha*CUB(coordinate[i]); 


else { 

acceleration = gamma|[i]*(coordinate[i+ 1]+coordinate[i-1] 
-2*coordinate[i])-beta*momentum[I] 

-(SQR(omega0[i]) + 2*eta*cos(2*omega*model_time-delta))*coordinate[1] 
+ alpha*CUB(coordinate[1]); 


old_momentum[i]=momentum{]]; 
momentum[1i]=momentum[i] + acceleration*time_int; 
\ /* DOFOR */ 

DOFOR(i,no_pendulums) { 
oldold_coordinate[i]=old_coordinate[i]; 
old_coordinate[1]=coordinate[1]; 
if(pinned_elements[i] = =0) { 
coordinate[i] = coordinate[i] +momentum[i]*time_ int; 


} 
3. Lower Cutoff, End Driven, Equation of Motion 


eqn motion() { 
/* The following loop calculates a round of velocities for one 
iteration of the model clock. This is where the equations 
of motion need to be located, if one wishes to change the 
model’s physics! */ 
int i,j,k; 
DOFOR(i,no_pendulums) 
if(coordinate[i] > 1000) { 
retrieve State(); 


DOFOR(i,no_pendulums) { 
if(i= =O) { 
acceleration = gamma[0]*(coordinate[ 1 ]-coordinate[0]) 
-beta*momentum[0] 
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-((coordinate[i]/sqrt(1_ + 2*SQR(coordinate[i]))) +eta*cos(omega*model_time)); 


else if(i= =no_pendulums-1){ 

acceleration = gamma|i]*(coordinate[i- 1 ]-coordinate[i]) 
-beta*momentum[!] 

-(coordinate[i]/sqrt(1 + 2*SQR(coordinate[i]))); 


else { 

acceleration = gamma|i]*(coordinate[i+ 1] +coordinate[i-1] 
-2*coordinate[i])-beta*momentum[!] 

-(coordinate[i]/sqrt(1 + 2*SQR(coordinate[i]))); 


old_momentum[i]=momentun 1]; 
momentum[i]=momentum|[i]+acceleration*time_int; 
} /* DOFOR */ 

DOFOR(i,no_pendulums) { 
old_coordinate[i] = coordinate[i]; 
if(pinned_elements[i]= =0) { 
coordinate[i] =coordinate[i] + momentum[i]*time_int; 


} 


4. Upper Cutoff, End Driven, Equation of Motion 


eqn motion( { 
/* The following loop calculates a round of velocities for one 
iteration of the model clock. This is where the equations 
of motion need to be located, if one wishes to change the 
model’s physics! */ 
int i,),k; 
DOFOR(i,no_pendulums) 
if(coordinate[i] > 100) { 
retrieve state(Q); 


DOFOR(i,no_pendulums) { 
if(i= =O) { 
acceleration = gamma[0]*(coordinate[ 1 ]-coordinate(0}) 
-beta*momentum[0] 
-(SQR(omega0[0])*coordinate[0] + eta*cos(omega*model_time))+alpha*CUB(coordinate[0)); 


else if(i= =no_pendulums-1){ 

acceleration = gamma|i]*(coordinate[i-1]-coordinate[i]) 
-beta*momentum[I] 

-(SQR(omega0[i]) *coordinate[i]) + alpha*CUB(coordinate[i]); 
} 


else { 
acceleration = gamma|i]*(coordinate[i+ 1]+coordinate[i-1] 
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-2*coordinate[1i])-beta*momentum[1] 
-(SQR(omega0[i])*coordinate[i]) + alpha*CUB(coordinate[i]}); 


old_momentum[i]=momentun{i]; 
momentum[i]=momentum[i] + acceleration*time_int; 
\ /* DOFOR */ 

DOFOR(i,no_pendulums) { 
old_coordinate[i] = coordinate[i]; 
if(pinned_elements[i]==0) { 
coordinate[i] =coordinate[i] +momentum[i]*time_int; 


} 
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